Para replicar las secciones [1.] OpenStreetMap y [2.] Otras fuentes de datos, debe descargar el siguiente proyecto de R y abrir el archivo clase-03.Rproj.

[1.] OpenStreetMap

OpenStreetMap (OSM) es un proyecto de mapeo de acceso abierto global, gratuito y licenciado bajo ODbL Licence. OSM recopila información geográfica capturada con dispositivos GPS móviles, ortofotografías y otras fuentes de información libre.

1.0 Configuración inicial

## Llamar pacman (contiene la función p_load)
require(pacman) 

## Llama/instala-llama las librerías listadas
p_load(tidyverse,rio,skimr,
       viridis, # paletas de colores
       sf, # Leer/escribir/manipular datos espaciales
       leaflet, # Visualizaciones dinámicas
       tmaptools, # geocode_OSM()
       osmdata) # Get OSM's data

1.1. Geocodificar direcciones

La función geocode_OSM() de la librería tmaptools se conecta a la API de OpenStreetMap y retorna el vértice con la coordenada geográfica del sitio/dirección buscado.

## Buscar un lugar público por el nombre
geocode_OSM("Casa de Nariño, Bogotá")
## $query
## [1] "Casa de Nariño, Bogotá"
## 
## $coords
##          x          y 
## -74.077549   4.595506 
## 
## $bbox
##       xmin       ymin       xmax       ymax 
## -74.077962   4.595103 -74.077078   4.595872

Puede adicionar el argumento as.sf=T para generar un objeto de clase sf:

## geocode_OSM no reconoce el caracter #, en su lugar se usa %23% 
cbd <- geocode_OSM("Centro Internacional, Bogotá", as.sf=T) 
cbd
## Simple feature collection with 1 feature and 7 fields
## Active geometry column: point
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -74.07057 ymin: 4.613615 xmax: -74.07057 ymax: 4.613615
## Geodetic CRS:  WGS 84
##                          query      lat       lon  lat_min  lat_max   lon_min
## 1 Centro Internacional, Bogotá 4.613615 -74.07057 4.613373 4.613889 -74.07098
##     lon_max                           bbox                      point
## 1 -74.07023 POLYGON ((-74.07098 4.61337... POINT (-74.07057 4.613615)

Puede visualizar el objeto point usando la librería leaflet:

## la función addTiles adiciona la capa de OpenStreetMap
leaflet() %>% addTiles() %>% addCircles(data=cbd)

1.2. Librería osmdata

osmdata es una librería de R que permite descargar y usar datos de OpenStreetMap (OSM).

1.2.1. Features disponibles

Puede acceder a la lista de features disponibles en OSM aquí. En R puede obtener un vector con los nombres de los features usando la función available_features():

available_features() %>% head(20)
##  [1] "4wd_only"                "abandoned"              
##  [3] "abutters"                "access"                 
##  [5] "addr"                    "addr:city"              
##  [7] "addr:conscriptionnumber" "addr:country"           
##  [9] "addr:county"             "addr:district"          
## [11] "addr:flats"              "addr:full"              
## [13] "addr:hamlet"             "addr:housename"         
## [15] "addr:housenumber"        "addr:inclusion"         
## [17] "addr:interpolation"      "addr:place"             
## [19] "addr:postbox"            "addr:postcode"

Cada feature contiene una lista de tags. Puede acceder al vector de tags usando la función available_tags(). Por ejemplo, puede acceder a la lista de amenity así:

available_tags("amenity") %>% head(20)
##  [1] "animal_boarding"        "animal_breeding"        "animal_shelter"        
##  [4] "animal_training"        "arts_centre"            "atm"                   
##  [7] "baby_hatch"             "baking_oven"            "bank"                  
## [10] "bar"                    "bbq"                    "bench"                 
## [13] "bicycle_parking"        "bicycle_rental"         "bicycle_repair_station"
## [16] "biergarten"             "boat_rental"            "boat_sharing"          
## [19] "brothel"                "bureau_de_change"

1.3. Descargar features

Para descargar features desde OSM, primero debe definir un espacio geográfico y obtener la caja de de coordenadas que lo contiene:

## obtener la caja de coordenada que contiene el polígono de Bogotá
opq(bbox = getbb("Bogotá Colombia"))
## $bbox
## [1] "4.4711754,-74.2235137,4.8331695,-74.0102577"
## 
## $prefix
## [1] "[out:xml][timeout:25];\n(\n"
## 
## $suffix
## [1] ");\n(._;>;);\nout body;"
## 
## $features
## NULL
## 
## $osm_types
## [1] "node"     "way"      "relation"
## 
## attr(,"class")
## [1] "list"           "overpass_query"
## attr(,"nodes_only")
## [1] FALSE

Puede almacenar un objeto osm de clase list y overpass_query con las estaciones de autobús para la ciudad de Bogotá:

## objeto osm
osm = opq(bbox = getbb("Bogotá Colombia")) %>%
      add_osm_feature(key="amenity" , value="bus_station") 
class(osm)
## [1] "list"           "overpass_query"

Puede acceder a los Simple Features del objeto osm usando la función osmdata_sf(). El objeto contiene una lista de objetos con los puntos, líneas y polígonos disponibles.

## extraer Simple Features Collection
osm_sf = osm %>% osmdata_sf()
osm_sf
## Object of class 'osmdata' with:
##                  $bbox : 4.4711754,-74.2235137,4.8331695,-74.0102577
##         $overpass_call : The call submitted to the overpass API
##                  $meta : metadata including timestamp and version numbers
##            $osm_points : 'sf' Simple Features Collection with 2708 points
##             $osm_lines : 'sf' Simple Features Collection with 102 linestrings
##          $osm_polygons : 'sf' Simple Features Collection with 381 polygons
##        $osm_multilines : NULL
##     $osm_multipolygons : 'sf' Simple Features Collection with 102 multipolygons

Puede acceder a los elementos de la lista y crear un objeto de clase sf:

## Obtener un objeto sf
bus_station = osm_sf$osm_points %>% select(osm_id,amenity) 
bus_station
## Simple feature collection with 2708 features and 2 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -74.22302 ymin: 4.530763 xmax: -74.03537 ymax: 4.812363
## Geodetic CRS:  WGS 84
## First 10 features:
##              osm_id amenity                   geometry
## 261630505 261630505    <NA> POINT (-74.14692 4.631477)
## 261630506 261630506    <NA> POINT (-74.14694 4.631353)
## 261630508 261630508    <NA> POINT (-74.14423 4.631102)
## 261630509 261630509    <NA> POINT (-74.14422 4.631225)
## 266537186 266537186    <NA> POINT (-74.11613 4.654154)
## 266537187 266537187    <NA> POINT (-74.11531 4.652883)
## 266537189 266537189    <NA> POINT (-74.11439 4.652935)
## 266537456 266537456    <NA>  POINT (-74.11531 4.65398)
## 266537457 266537457    <NA>  POINT (-74.11507 4.65568)
## 292774696 292774696    <NA> POINT (-74.14562 4.631682)

Visualiza el objeto bus_station:

## Pintar las estaciones de autobus
leaflet() %>% addTiles() %>% addCircleMarkers(data=bus_station , col="red")

[2.] Otras fuentes de datos

2.1 Censo Nacional

Aquí puede acceder a los datos del Censo Nacional de Población de 2018 para Colombia.

## censo data
browseURL("https://microdatos.dane.gov.co//catalog/643/get_microdata")

##=== variables ===##

## COD_DANE_ANM: Codigo DANE de manzana
## UA_CLASE: ID
## COD_ENCUESTAS: ID de encuesta
## U_VIVIENDA: ID de vivienda
## H_NRO_CUARTOS: Número de cuartos en total
## HA_TOT_PER: Total personas en el hogar
## V_TOT_HOG: Total de hogares en la vivienda
## VA1_ESTRATO: Estrato de la vivienda (según servicio de energía)

##=== load data ===##

## unzip file
unzip(zipfile="input/11_BOGOTA_CSV.zip" , exdir="input/." , overwrite=T) 

## data manzanas
mgn <- import("input/censo_2018/CNPV2018_MGN_A2_11.CSV")
colnames(mgn)
mgn <- mgn %>% select(COD_DANE_ANM,UA_CLASE,COD_ENCUESTAS,U_VIVIENDA)

## data hogar
hog <- import("input/censo_2018/CNPV2018_2HOG_A2_11.CSV")
colnames(hog)
hog <- hog %>% select(UA_CLASE,COD_ENCUESTAS,U_VIVIENDA,H_NROHOG,H_NRO_CUARTOS,HA_TOT_PER)

## data vivienda
viv <- import("input/censo_2018/CNPV2018_1VIV_A2_11.CSV") 
colnames(viv)
viv <- viv %>% select(COD_ENCUESTAS,UA_CLASE,U_VIVIENDA,V_TOT_HOG,VA1_ESTRATO)

## join hogar-vivienda
viv_hog <- left_join(hog,viv,by=c("COD_ENCUESTAS","U_VIVIENDA","UA_CLASE"))

## joing mnz-hogar-vivienda
viv_hog_mgn <- left_join(viv_hog,mgn,by=c("UA_CLASE","COD_ENCUESTAS","U_VIVIENDA"))

##=== collapse data ===##
db <- viv_hog_mgn %>%
      group_by(COD_DANE_ANM) %>% 
      summarise(med_H_NRO_CUARTOS=median(H_NRO_CUARTOS,na.rm=T), 
                sum_HA_TOT_PER=sum(HA_TOT_PER,na.rm=T), 
                med_V_TOT_HOG=median(V_TOT_HOG,na.rm=T),
                med_VA1_ESTRATO=median(VA1_ESTRATO,na.rm=T))

## export data
export(db,"output/mnz_censo_2018.rds")

## delete files
unlink("input/censo_2018/",recursive = T)
unlink("input/__MACOSX/",recursive = T)

2.2 MGN

Aquí puede acceder al Marco Geoestadístico Nacional para Colombia.

## MGN data
browseURL("https://geoportal.dane.gov.co/servicios/descarga-y-metadatos/descarga-mgn-marco-geoestadistico-nacional/")

## load data
mgn <- st_read("input/mgn/MGN_URB_MANZANA.shp")
## Reading layer `MGN_URB_MANZANA' from data source 
##   `/Users/eduard-martinez/Dropbox/teaching/bid-data-real-state/2023-01/clase-03/public/input/mgn/MGN_URB_MANZANA.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 43438 features and 13 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -74.36463 ymin: 3.888847 xmax: -74.00074 ymax: 4.833006
## Geodetic CRS:  WGS 84
mgn %>% head()
## Simple feature collection with 6 features and 13 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -74.13846 ymin: 4.661075 xmax: -74.11615 ymax: 4.683528
## Geodetic CRS:  WGS 84
##   DPTO_CCDGO MPIO_CCDGO CLAS_CCDGO SETR_CCDGO SECR_CCDGO CPOB_CCDGO SETU_CCDGO
## 1         11      11001          1        000         00   11001000       6402
## 2         11      11001          1        000         00   11001000       6402
## 3         11      11001          1        000         00   11001000       6403
## 4         11      11001          1        000         00   11001000       6403
## 5         11      11001          1        000         00   11001000       6302
## 6         11      11001          1        000         00   11001000       6302
##   SECU_CCDGO MANZ_CCDGO             MANZ_CCNCT MANZ_CAG  Shape_Leng
## 1         03         23 1100110000000064020323   145420 0.006490357
## 2         09         19 1100110000000064020919   145102 0.004081371
## 3         02         15 1100110000000064030215   144714 0.002126900
## 4         02         17 1100110000000064030217   144665 0.001641266
## 5         04         10 1100110000000063020410   147399 0.002941721
## 6         02         03 1100110000000063020203   147590 0.003051028
##     Shape_Area                       geometry
## 1 1.413612e-06 POLYGON ((-74.1378 4.677946...
## 2 9.998505e-07 POLYGON ((-74.13436 4.67957...
## 3 2.572230e-07 POLYGON ((-74.13465 4.68263...
## 4 1.500361e-07 POLYGON ((-74.1348 4.683005...
## 5 4.549105e-07 POLYGON ((-74.1215 4.663893...
## 6 6.541770e-07 POLYGON ((-74.11667 4.66107...
table(mgn$CLAS_CCDGO)
## 
##     1     2 
## 43361    77
mgn <- mgn %>% subset(CLAS_CCDGO==1)
mgn <- mgn %>% select(MANZ_CCNCT)

2.3 CENSO-MGN

## load data censo 2018
censo <- import("output/mnz_censo_2018.rds")
censo
## # A tibble: 38,567 × 5
##    COD_DANE_ANM           med_H_NRO_CUARTOS sum_HA_TOT_PER med_V_TOT_HOG med_V…¹
##    <chr>                              <dbl>          <int>         <dbl>   <dbl>
##  1 1100110000000000000000                NA              0            NA      NA
##  2 1100110000000011010101                 6             13             1       2
##  3 1100110000000011010102                 3            442             1       2
##  4 1100110000000011010103                 4            787             1       2
##  5 1100110000000011010104                 4           2511             1       2
##  6 1100110000000011010105                 2            224             1       2
##  7 1100110000000011010106                 3            223             1       2
##  8 1100110000000011010107                 3            279             1       2
##  9 1100110000000011010108                 3            121             1       2
## 10 1100110000000011020101                 4            166             1       2
## # … with 38,557 more rows, and abbreviated variable name ¹​med_VA1_ESTRATO
## agregar a manzanas de mgn info de censo
mnz_censo <- left_join(mgn,censo,by=c("MANZ_CCNCT"="COD_DANE_ANM"))
mnz_censo
## Simple feature collection with 43361 features and 5 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -74.22353 ymin: 4.46292 xmax: -74.00074 ymax: 4.833006
## Geodetic CRS:  WGS 84
## First 10 features:
##                MANZ_CCNCT med_H_NRO_CUARTOS sum_HA_TOT_PER med_V_TOT_HOG
## 1  1100110000000064020323                 4            465             1
## 2  1100110000000064020919                 4            656             1
## 3  1100110000000064030215                NA             NA            NA
## 4  1100110000000064030217                 4             93             1
## 5  1100110000000063020410                 6             55             1
## 6  1100110000000063020203                 4            538             1
## 7  1100110000000063020239                 5            175             1
## 8  1100110000000063020517                 4            647             1
## 9  1100110000000063020513                 4            439             1
## 10 1100110000000063120117                 5            104             1
##    med_VA1_ESTRATO                       geometry
## 1                3 POLYGON ((-74.1378 4.677946...
## 2                3 POLYGON ((-74.13436 4.67957...
## 3               NA POLYGON ((-74.13465 4.68263...
## 4                3 POLYGON ((-74.1348 4.683005...
## 5                4 POLYGON ((-74.1215 4.663893...
## 6                4 POLYGON ((-74.11667 4.66107...
## 7                4 POLYGON ((-74.11826 4.66223...
## 8                4 POLYGON ((-74.11716 4.66048...
## 9                4 POLYGON ((-74.11905 4.66280...
## 10               4 POLYGON ((-74.1238 4.668334...
## visualizar estratos
ggplot() + geom_sf(data=mnz_censo , aes(fill=as.factor(med_VA1_ESTRATO)) , col=NA)

## visualizar densidad poblacional
ggplot() + geom_sf(data=mnz_censo , aes(fill=sum_HA_TOT_PER) , col=NA) +
scale_fill_viridis(option="inferno")

## exportar sf de manzanas con info del censo
export(mnz_censo,"output/mgn_censo_2018.rds")

Referencias

  • Lovelace, R., Nowosad, J., & Muenchow, J. (2019). Geocomputation with R. [Ver aquí]

    • Cap. 4: Spatial Data Operations
    • Cap. 5: Geometry Operations
    • Cap. 6: Reprojecting geographic data
    • Cap. 11: Statistical learning

  • Bivand, R. S., Pebesma, E. J., Gómez-Rubio, V., & Pebesma, E. J. (2013). Applied spatial data analysis with R. [Ver aquí]

    • Cap. 4: Spatial Data Import and Export
    • Cap. 7: Spatial Point Pattern Analysis
    • Cap. 8: Interpolation and Geostatistics